BEHAVIORAL DECISION THEORY (BDT) EXPERIMENTAL DESIGN, MODEL ESTIMATION AND SELECTION

Data Science Conference 2022 Tech Tutorial, 15. November 2022.

Abstract

We introduce the fundamental concepts, typical approaches to experimental design and measurement, and parametric model estimation via MLE in Behavioral Decision Theory (BDT). In the following 1.5 hours of this tech tutorial in the scope of DSC 2022 we will consider only choice under risk (and not uncertainty), and study two important models. We will review and provide tools in R for the MLE estimation of the (normative) Expected Utility Theory (EUT, von Neumann–Morgenstern) and the (descriptive) Cumulative Prospect Theory (Tversky-Kahneman). We approach the model selection in BDT via k-fold cross-validation, more characteristic of the ML community and rarely used in empirical approaches in microeconomics, behavioral economics, and cognitive psychology. However, typical, ordinary approaches to k-fold CV fail because of the deep dependence between the experimental design (i.e. the structure of the data set) and model estimation procedures. Understanding the fundaments of BDT and being able to estimate parametric models of choice under risk is a prerequisite for the design of any behavioral intervention, an aproach that gained so much popularity and attracted so much attention following the recognition of the central importance of behavioral economics (e.g Kahneman’s (2002) and Thaler’s (2017) Nobel Memorial Prize in Economic Sciences) to the study of human decision making. Beyond, characterizations of individual decision makers by estimated BTD model parameters provide new sets of features for any further risk modeling.

Goran S. Milovanović, Phd DataKolektiv, Owner Smartocto, Senior Data Scientist

---

Feedback should be send to goran.milovanovic@datakolektiv.com. These notebook accompanies the BEHAVIORAL DECISION THEORY (BDT) EXPERIMENTAL DESIGN, MODEL ESTIMATION AND SELECTION Data Science Conference 2022 Tech Tutorial, 15. November 2022, 16:30 - 18:00 CET

Setup

# - setup
data_dir <- paste0(getwd(), "/_data/")
meq_data_dir <- paste0(data_dir, "_meq/")
kt92_data_dir <- paste0(data_dir, "_kt92/")

# - libs
library(tidyverse)
library(scales)
library(plotly)
library(snowfall)

Fundamental concepts in Behavioral Decision Theory

Risk aversion and risk attitudes

The risky prospect and the sure thing have the same expected value:

\[E[Risky Prospect] = E[X=100, p=.5;X=0, p=.5] = .5\times100 + .5\times0 = 50\]

\[E[Sure Thing] = 50\]

Why are people risk averse?

Utility Functions

Utility Function of an individual decision maker in choice under risk

Daniel Bernoulli, “Specimen Theoriae Novae de Mensura Sortis”, Comentarii Academiae Scientiarum Imperialis Petroolitanae, Tomus V, 1738, pp. 175-192

Risk aversion follows as a consequences of the characteristics of the decision maker’s utility function.

The power utility function:

\[u(x) = x^\rho\] The decision maker is risk averse for \(0 < \rho < 1\) and risk seeking for \(\rho>1\); risk seeking is not seen as a rational risk attitude.

# - power utility function over gains
x <- seq(1, 100, .1)
rho = .67
ux <- x^rho
ra_frame <- data.frame(value = x,
                       utility = ux)
ggplot(data = ra_frame,
       aes(x = value, y = utility)) + 
  geom_point(size = .25, color = "blue") +
  geom_segment(aes(x = 50, 
                   y = 0, 
                   xend = 50, 
                   yend = 50^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = 50^rho, 
                   xend = 50, 
                   yend = 50^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 100, 
                   y = 0, 
                   xend = 100, 
                   yend = 100^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = 100^rho, 
                   xend = 100, 
                   yend = 100^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = 100^rho, 
                   xend = 100, 
                   yend = 100^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = (100^rho)/2, 
                   xend = 100, 
                   yend = (100^rho)/2),
               colour = "red",
               linetype = "dashed",
               size = .1) + 
  ggtitle("Power utility function, rho=.67, Risk Aversion") +
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

# - power utility function over gains
x <- seq(1, 100, .1)
rho = 1.2
ux <- x^rho
ra_frame <- data.frame(value = x,
                       utility = ux)
ggplot(data = ra_frame,
       aes(x = value, y = utility)) + 
  geom_point(size = .25, color = "red") +
  geom_segment(aes(x = 50, 
                   y = 0, 
                   xend = 50, 
                   yend = 50^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = 50^rho, 
                   xend = 50, 
                   yend = 50^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 100, 
                   y = 0, 
                   xend = 100, 
                   yend = 100^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = 100^rho, 
                   xend = 100, 
                   yend = 100^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = 100^rho, 
                   xend = 100, 
                   yend = 100^rho),
               colour = "black",
               linetype = "dashed",
               size = .1) + 
  geom_segment(aes(x = 0, 
                   y = (100^rho)/2, 
                   xend = 100, 
                   yend = (100^rho)/2),
               colour = "red",
               linetype = "dashed",
               size = .1) + 
  ggtitle("Power utility function, rho=1.2, Risk Seeking") +
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

Daniel Bernoulli, 1738. People make choices under risk by following the principle of the Maximum Expected Utility (MEU):

\[EU(Prospect) = \sum_{i=1}^{N}p_iu(x_i)\]

The approach was axiomatized in 1947. only by John von Neumann and Oskar Morgenstern in their “Theory of Games and Economic Behavior”, by providing a set of axioms that introduce the preference relation and a proof of a representation theorem showing that any decision maker who obeys the axioms of rational choice is at the same time a EU maximizer (and vice versa). Since then: the von Neumann-Morgenstern utility (or vNM Expected Utility Theory).

Kahneman and Tversky’s Reference-Dependent Utility Function

People are risk seeking the domain of losses.

Example. The “Asian disease” framing problem (Tversky & Kahneman, 1981):

Imagine that the US is preparing for the outbreak of an unusual Asian disease, which is expected to kill 600 people. Two alternative programs to combat the disease have been proposed.

• Program A: 200 people will be saved (72% of participants expressed a preference for A).
• Program B: 1/3 probability that 600 people will be saved, and a 2/3 probability that no people will be saved.

This finding illustrates risk aversion in the domain of gains.

• Program C: 400 people will die (22% opted for C).
• Program D: 1/3 probability that nobody will die, and a 2/3 probability that 600 people will die.

This finding illustrates risk seeking in the domain of losses.

The Reference-Dependent Utility Function:

# - power reference-dependent utility function over gains and losses
x <- seq(-100, 100, .1)
rho = .67
ux <- ifelse(x>=0, x^rho, -(-x)^rho)
ra_frame <- data.frame(value = x,
                       utility = ux)
ggplot(data = ra_frame,
       aes(x = value, y = utility)) + 
  geom_point(size = .25, color = "blue") + 
  geom_hline(yintercept=0, size = .1, colour = "black", linetype = "dashed") +
  geom_vline(xintercept=0, size = .1, colour = "black", linetype = "dashed") +
  ggtitle("Reference-Dependent Power Utility function, rho=.67") +
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

Loss aversion

Kahneman, Knetch & Thaler experiment, 1991.

Half of the participants in the experimental study receive something for free at the beginning of the session, e.g. a beatiful cup

For those who own the item: place the offer.

$7 and 12 cents

• For those who do not own the item: place the bid.

2$ and 87 cents

That would be a 2.48 ratio in favor of the offered priced.

We need another correction of the utility function. Enter Loss Aversion, \(\lambda\):

# - power reference-dependent utility function over gains and losses
x <- seq(-100, 100, .1)
rho = .67
lambda = 2.2
ux <- ifelse(x>=0, x^rho, -lambda*(-x)^rho)
ra_frame <- data.frame(value = x,
                       utility = ux)
ggplot(data = ra_frame,
       aes(x = value, y = utility)) + 
  geom_point(size = .25, color = "blue") + 
  geom_hline(yintercept=0, size = .1, colour = "black", linetype = "dashed") +
  geom_vline(xintercept=0, size = .1, colour = "black", linetype = "dashed") +
  ggtitle("Loss Averse Reference-Dependent Power Utility function, \nrho=.67, lambda=2.2") +
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

If \(\lambda > 1\), we observe loss aversion; however, if \(0<\lambda<1\), we observe gain seeking, and that happens in empirical settings more often than people think.

But there is more to it.

Probability weighting

The four-fold pattern of risk attitudes (Kahneman & Tversky), e.g.

• If a risky prospect offers only a 5% chance for EUR 100 and a 95% chance to win nothing – people tend to be risk seeking (i.e. asking for more than EUR 5 in exchange for the lottery), but
• If a risky prospect offers a 95% chance for EUR 100 and only a 5% chance to win nothing – people tend to be risk averse (i.e. asking less than EUR 95in exchange for the lottery);

and of course we observe the reflection effect in the domain of losses:

• If a risky prospect offers only a 5% chance for a loss of EUR 100 and a 95% chance to lose nothing – people tend to be risk averse (i.e. offering more than EUR 5 to avoid the risky event), but
• If a risky prospect offers a 95% chance to lose EUR 100 and only a 5% chance to lose nothing – people tend to be risk seeking (i.e. offering less than EUR 95 to avoid the risky event).

It is like if there was a probability weighting function, w(p) such that \(w(p) > p\) for small \(p\) and \(w(p) < p\) for high \(p\), with an inflection point somewhere. Many functional forms were proposed for \(w(p)\); we will mention only the Prelec’s one-paramater form (1988):

\[w(p)=exp(-(-ln(p))^\gamma)), 0<\gamma<1 \]

# - Prelec's single-parameter probability weighting function
p_weight <- function(p, gamma) {
  wp <- exp(-(-log(p))^gamma)
  return(wp)
} 
p <- seq(0,1,.01)
gamma <- .65
wp <- p_weight(p, gamma)
wp_frame <- data.frame(probability = p,
                       w_p = wp)
ggplot(data = wp_frame,
       aes(x = probability, y = wp)) + 
  geom_path(size = .25, color = "blue", group=1) + 
  ggtitle("Prelec's one-parameter w(p), gamma=.65") +
  geom_segment(aes(x = 0, 
                   y = 0, 
                   xend = 1, 
                   yend = 1),
               colour = "black",
               linetype = "dashed",
               size = .1) +
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

Prospect Theory

Finally,

\[u_g(x)=x^{\rho_g}\]

\[u_l(x)=-{\lambda}x^{\rho_l}\]

\[w_g(p)=exp(-(-ln(p_g))^{\gamma_g}))\] \[w_l(p)=exp(-(-ln(p_l))^{\gamma_l}))\]

And then in the simplest case we have

\[PT_G(Prospect) = \sum_{i=1}^{N}w_g(p_i)u_g(x_i)\]

for gains and

\[PT_L(Prospect) = \sum_{i=1}^{N}w_l(p_i)u_l(x_i)\] for losses. For mixed prospects (i.e. lotteries including both gains and losses), we sum up the \(PT()\) for the positive and the negative part:

\[PT_{MIXED}(Prospect) = PT_G(Prospect)+PT_L(Prospect)\]

Measurement and Data

Inspect some Monetary Equivalents (MEq) data (Kahneman & Tversky, 1992):

# - load some MEq experimental data
meq_data <- read.csv(paste0(kt92_data_dir, "kt92.csv"),
                     header = TRUE,
                     check.names = FALSE,
                     row.names = 1,
                     stringsAsFactors = FALSE)
# - percents to probability in meq_data
meq_data$p1 <- meq_data$p1/100
meq_data$p2 <- meq_data$p2/100
print(meq_data)

Estimating the Expected Utility Theory (EUR) model via LS

Start simple: let’s estimate the EUT parameter (there is only one: the power-utility \(\rho\) exponent) via Least Squares.

Utility function

utility_power() is our EUT utility function.

# - power utility function
utility_power <- function(x, rho) {
  u <- x^rho
  return(u)
}

Prediction

A function eut_predict() to predict monetary equivalents from a dataset via EUT:

# - function: EUT predictions
eut_predict <- function(data, rho) {
  
  # - utility functions
  u1 <- ifelse(data$v1 >= 0,
               utility_power(data$v1, rho),
               -utility_power(-data$v1, rho)
               )
  u2 <- ifelse(data$v2 >= 0,
               utility_power(data$v2, rho),
               -utility_power(-data$v2, rho)
               )
  
  # - predictions: utility scale
  predictions <- data$p1*u1 + data$p2*u2
  
  # - predictions: value scale
  predictions <- ifelse(predictions >= 0,
                        predictions^(1/rho),
                        -((-predictions)^(1/rho)))
  
  # - output
  return(predictions)
  
}

NOTE. The predictions: value scale part of the code is very important, because we do not wish to make predictions on the utility scale but rather on the value (MEq) scale. So, let’s not forget that

\[MEq = u(x)^{1/\rho}\] for gains, and

\[MEq = -(u(-x)^{1/\rho})\] for losses.

And test our predictions for \(\rho=.76\):

preds <- eut_predict(data = meq_data, rho = .76)
print(preds)

 [1]    2.4164651   20.0852807   43.5274054   -2.4164651  -20.0852807  -43.5274054
 [7]    1.9414152   16.1367400   40.1705613   68.4868080   93.4735968   -1.9414152
[13]  -16.1367400  -40.1705613  -68.4868080  -93.4735968    0.4671443    9.6658605
[19]   80.3411227  174.1096214  197.3725857   -0.4671443   -9.6658605  -80.3411227
[25] -174.1096214 -197.3725857    0.9342886  394.7451714   -0.9342886 -394.7451714
[31]   54.6116284   73.9745732   94.6469643  -54.6116284  -73.9745732  -94.6469643
[37]   54.3353025   72.4912535   96.8153714  122.7129797  144.4385997  -54.3353025
[43]  -72.4912535  -96.8153714 -122.7129797 -144.4385997  104.5872226  123.4136252
[49]  147.9491464  173.5050929  194.6292955 -104.5872226 -123.4136252 -147.9491464
[55] -173.5050929 -194.6292955

SSE: objective

Now the objective function eut_sse(): it computes the \(SSE\) following the predictive pass of eut_predict()

# - optimize the EUT model via LS
eut_sse <- function(params) {
  rho <- abs(params[1])
  preds <- eut_predict(meq_data, rho)
  sse = sum((meq_data$meq - preds)^2)
  return(sse)
}

Test for \(\rho=.76\):

eut_sse(params=.76)

[1] 7535.473

EUT model optimization

And we want to minimize eut_sse() for the MEq dataset at hand. We will use Nelder-Mead with R optim():

# - random initial value
init_rho <- runif(1, 0, 1)
solution <- optim(par = init_rho, 
                  fn = eut_sse, 
                  method = "Nelder-Mead", 
                  control = list("maxit" = 10e6))
print(paste0("Estimated Rho of ", round(solution$par, 2)))

[1] "Estimated Rho of 0.7"

print(paste0("Convergence check: ", solution$convergence))

[1] "Convergence check: 0"

Now let’s predict the observed MEqs from the estimated parameter:

# - predict from the optimized EUT model
meq_data$eut_meq <- eut_predict(meq_data,
                                rho = solution$par)
print(meq_data)

And plot the predictions:

# - plot predictions
ggplot(data = meq_data,
       aes(x = eut_meq, y = meq)) + 
  geom_smooth(method = "lm", size = .25, color = "red", se = FALSE) + 
  geom_point(size = 1.5, color = "black") + 
  geom_point(size = 1, color = "white") + 
  ggtitle(paste0("EUT, Estimated Rho of ", round(solution$par, 2))) +
  xlab("Predicted MEq") + 
  ylab("Observed MEq") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

Let’s run 100 optimizations to check if we are stuck in some local minimum only:

# - many optimization runs
init_rho <- runif(100, 0, 1)

# - run
results <- lapply(init_rho, function(x) {
  
  # - optimization run
  optim_sol <- tryCatch(optim(par = x,
                              fn = eut_sse,
                              method = "Nelder-Mead",
                              control = list("maxit" = 10e6)),
                        error = function(condition) {
                          return(NULL)
                        })
  
  if (!is.null(optim_sol)) {
    
    # - predictions
    preds <- eut_predict(meq_data, rho = optim_sol$par)
    
    return(
      list(rho = optim_sol$par,
           sse = optim_sol$value,
           convergence = optim_sol$convergence,
           predictions = preds)
      )
    
  } else {
      
      return(
      list(rho = NULL,
           sse = NULL,
           convergence = NULL,
           predictions = NULL)
      )
    
  }
  
  
})

# - find best solution
sse_list <- sapply(results, function(x) x$sse)
best_run <- which.min(sse_list)[1]
print(paste0("Check convergence: ", results[[best_run]]$convergence))

[1] "Check convergence: 0"

estimated_rho <- results[[best_run]]$rho
optimal_predictions <- results[[best_run]]$predictions

# - plot predictions
plot_frame <- meq_data
meq_data$eut_meq <- optimal_predictions
ggplot(data = meq_data,
       aes(x = meq, y = eut_meq)) + 
  geom_smooth(method = "lm", size = .25, color = "red", se = FALSE) + 
  geom_point(size = 1.5, color = "black") + 
  geom_point(size = 1, color = "white") + 
  ggtitle(paste0("EUT, Estimated Rho of ", round(estimated_rho, 2))) +
  xlab("Predicted MEq") + 
  ylab("Observed MEq") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

EUT model error surface

The EUT model error surface from a sample of 100,000 random parameter values:

# - sample parameters from a range for plotting purposes
sample_parameters <- data.frame(rho = runif(100000, 0, 1))
sample_parameters$sse <- sapply(sample_parameters$rho, eut_sse)
head(sample_parameters)

ggplot(data = sample_parameters,
       aes(x = rho, y = sse)) + 
  geom_line(size = .25, color = "black", group=1) + 
  ggtitle("EUT Model Error Surface ") +
  xlab("Rho") + 
  ylab("SSE") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

Estimating the Expected Utility Theory (EUR) model via MLE

Approach to MLE for EUT

For the Maximum-Likelihood Estimation (MLE) of EUT, let’s assume the following:

• the decision maker computes the \(u(x)\) internally, and the
• responds with \(u(x)+Normal(\mu=0,\sigma)\), i.e. adding a random noise to their response.

Then we can obtain the likelihood of each MEq observation from \(Normal(response-prediction, \sigma)\), estimating one additional response model parameter \(\sigma\).

The function eut_likelihood returns the negative log-likelihood for a dataset given parameters \(\rho,\sigma\):

# - function: EUT model likelihood
eut_likelihood <- function(data, rho, sigma) {
  
  # - utility functions
  u1 <- ifelse(data$v1 >= 0,
               utility_power(data$v1, rho),
               -utility_power(-data$v1, rho)
               )
  u2 <- ifelse(data$v2 >= 0,
               utility_power(data$v2, rho),
               -utility_power(-data$v2, rho)
               )
  
  # - predictions: utility scale
  predictions <- data$p1*u1 + data$p2*u2
  
  # - predictions: value scale
  predictions <- ifelse(predictions >= 0,
                        predictions^(1/rho),
                        -((-predictions)^(1/rho)))
  
  # - likelihood
  # - assume that the decision maker computes 
  # - the utility of the gamble and its MEQ, and then
  # - responds with an added random noise of Normal(0, sigma)
  neg_loglike <- -sum(
    dnorm(predictions - data$meq, mean = 0, sd = sigma, log = TRUE)
  )

  # - output
  return(neg_loglike)
  
}

Test eut_likelihood for \(\rho=.76\) and \(\sigma=11\):

nll <- eut_likelihood(data = meq_data, rho = .76, sigma = 11)
print(nll)

[1] 216.881

The objective

Introduce the objective eut_mle() function

# - optimize the EUT model via LS
eut_mle <- function(params) {
  rho <- abs(params[1])
  sigma <- abs(params[2])
  nll <- eut_likelihood(data = meq_data, rho, sigma)
  return(nll)
}

Optimize EUT model

And optimize eut_mle() for meq_data:

# - random initial value
init_rho <- runif(1, 0, 1)
init_sigma <- runif(1, 0, 1)
solution <- optim(par = c(init_rho, init_sigma), 
                  fn = eut_mle, 
                  method = "Nelder-Mead", 
                  control = list("maxit" = 10e6))
print(paste0("Estimated Rho of ", round(abs(solution$par[1]), 2)))

[1] "Estimated Rho of 0.7"

print(paste0("Estimated Sigma of ", round(abs(solution$par[2]), 2)))

[1] "Estimated Sigma of 11.53"

print(paste0("Convergence check: ", solution$convergence))

[1] "Convergence check: 0"

# - predict from the optimized EUT model
meq_data$eut_meq <- eut_predict(meq_data,
                                rho = abs(solution$par[1]))

# - plot predictions
ggplot(data = meq_data,
       aes(x = meq, y = eut_meq)) + 
  geom_smooth(method = "lm", size = .25, color = "red", se = FALSE) + 
  geom_point(size = 1.5, color = "black") + 
  geom_point(size = 1, color = "white") + 
  ggtitle(paste0("EUT, Estimated Rho of ", round(abs(solution$par), 2))) +
  xlab("Predicted MEq") + 
  ylab("Observed MEq") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

Check by 100 optimization runs:

# - many optimization runs
init_pars <- lapply(1:100, function(x) {
  return(
    list(
      init_rho = runif(1, 0, 1),
      init_sigma = runif(1, 0, 1)
    )
  )
})

# - run
results <- lapply(init_pars, function(x) {
  
  init_rho <- x$init_rho
  init_sigma <- x$init_sigma
  
  # - optimization run
  optim_sol <- tryCatch(optim(par = c(init_rho, init_sigma),
                              fn = eut_mle,
                              method = "Nelder-Mead",
                              control = list("maxit" = 10e6)),
                        error = function(condition) {
                          return(NULL)
                        })
  
  if (!is.null(optim_sol)) {
    
    # - predictions
    preds <- eut_predict(meq_data, rho = abs(optim_sol$par[1]))

    return(
      list(rho = optim_sol$par[1],
           sigma = optim_sol$par[2],
           nll = optim_sol$value,
           convergence = optim_sol$convergence,
           predictions = preds)
      )
    
  } else {
      
      return(
      list(rho = NULL,
           sigma = NULL,
           nll = NULL,
           convergence = NULL,
           predictions = NULL)
      )
    
  }
  
  
})

# - find best solution
nll_list <- sapply(results, function(x) x$nll)
best_run <- which.min(nll_list)[1]
print(paste0("Check convergence: ", results[[best_run]]$convergence))

[1] "Check convergence: 0"

estimated_rho <- abs(results[[best_run]]$rho)
print(paste0("Estimated Rho: ", estimated_rho))

[1] "Estimated Rho: 0.699160388762866"

estimated_sigma <- abs(results[[best_run]]$sigma)
print(paste0("Estimated Sigma: ", estimated_sigma))

[1] "Estimated Sigma: 11.5171137364202"

optimal_predictions <- results[[best_run]]$predictions

# - plot predictions
plot_frame <- meq_data
meq_data$eut_meq <- optimal_predictions
ggplot(data = meq_data,
       aes(x = meq, y = eut_meq)) + 
  geom_smooth(method = "lm", size = .25, color = "red", se = FALSE) + 
  geom_point(size = 1.5, color = "black") + 
  geom_point(size = 1, color = "white") + 
  ggtitle(paste0("EUT, Estimated Rho of ", round(estimated_rho, 2))) +
  xlab("Predicted MEq") + 
  ylab("Observed MEq") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

The EUT Negative Log-Likelihood surface

# - sample parameters from a range for plotting purposes
sample_parameters <- data.frame(rho = runif(100000, 0, 1.5), 
                                sigma = runif(100000, 0, 15))
sample_parameters$nll <- apply(sample_parameters, 1, function(x) {
    eut_mle(params = c(x[1], x[2]))
})
head(sample_parameters)

plot_ly() %>% 
  add_trace(data = sample_parameters,  
            x = sample_parameters$rho, 
            y = sample_parameters$sigma, 
            z = sample_parameters$nll, 
            type = "mesh3d", 
            intensity = sample_parameters$nll) %>% 
  layout(
    modebar = list(orientation = "v"), 
    title = "The EUT Negative Log-Likehood function",
    scene = list(
      xaxis = list(title = "Rho", range = c(0, 1.5)),
      yaxis = list(title = "Sigma", range = c(0, 15)),
      zaxis = list(title = "Neg-LogLike", range = c(0, 1000))
    )) %>% hide_colorbar()

Estimating the (Cumulative) Prospect Theory (CPT) model via MLE

We will begin with a different, larger data set (Milovanovic, 2013):

# - load some MEq experimental data
meq_data <- read.csv(paste0(meq_data_dir, "es4v1.csv"),
                     header = TRUE,
                     check.names = FALSE,
                     row.names = 1,
                     stringsAsFactors = FALSE)
# - percents to probability in meq_data
meq_data$p1 <- meq_data$p1/100
meq_data$p2 <- meq_data$p2/100
print(meq_data)

The PT Reference-Dependent Power Utility Function

Enters rd_utility_power():

# - reference_dependent power utility function
rd_utility_power <- function(x, rho_g, rho_l, lambda) {
  u <- ifelse(x >= 0, 
              x^rho_g,
              -lambda*((-x)^rho_l))
  return(u)
}

Probability Weighting Function

Enters p_weight()

# - Prelec's single-parameter probability weighting function
# - NOTE: w(p) will be applied over gains and losses separately 
p_weight <- function(p, gamma) {
  
  wp <- exp(-(-log(p))^gamma)
  return(wp)

}

The Prospect Theory Negative Log-Likelihood

PT optimization and predictions

cpt_likelihood() returns the negative log-likelihood for a given dataset:

# - function: CPT model likelihood
cpt_likelihood <- function(data, 
                           rho_g, 
                           rho_l, 
                           lambda, 
                           gamma_gain, 
                           gamma_loss,
                           sigma) {
  
  # - utility functions
  u1 <- rd_utility_power(data$v1, rho_g, rho_l, lambda)
  u2 <- rd_utility_power(data$v2, rho_g, rho_l, lambda)
  
  # - probability weighting
  wp1 <- vector(mode = "numeric", length = length(data$v1))
  ix_gain_1 <- which(data$v1 >= 0)
  wp1[ix_gain_1] <- p_weight(data$p1[ix_gain_1], gamma_gain)
  ix_loss_1 <- which(data$v1 < 0)
  wp1[ix_loss_1] <- p_weight(data$p1[ix_loss_1], gamma_loss)
  wp2 <- 1-wp1
  
  # - predictions: utility scale
  predictions <- wp1*u1 + wp2*u2
  
  # - predictions: value scale
  predictions <- ifelse(predictions >= 0,
                        predictions^(1/rho_g),
                        -(((-predictions)/lambda)^(1/rho_l))
                        )
  
  # - likelihood
  # - assume that the decision maker computes 
  # - the utility of the gamble and its MEQ, and then
  # - responds with an added random noise of Normal(0, sigma)
  neg_loglike <- -sum(
    dnorm(predictions - data$meq, mean = 0, sd = sigma, log = TRUE)
  )

  # - output
  return(neg_loglike)
  
}

Test cpt_likelihood() for a set of parameters:

nll <- cpt_likelihood(data = meq_data, 
                      rho_g = .71, 
                      rho_l = .69, 
                      lambda = 2, 
                      gamma_gain = .67, 
                      gamma_loss = .61, 
                      sigma = .14)
print(nll)

[1] 5752202

And finally minimize cpt_likelihood() for meq_data:

# - optimize the CPT model via LS
cpt_mle <- function(params) {
  rho_g <- abs(params[1])
  rho_l <- abs(params[2])
  lambda <- abs(params[3])
  gamma_gain <- abs(params[4])
  gamma_loss <- abs(params[5])
  sigma <- abs(params[6])
  nll <- cpt_likelihood(data = meq_data, 
                        rho_g, rho_l, lambda,
                        gamma_gain, gamma_loss, sigma)
  return(nll)
}

# - random initial value
init_rho_g <- runif(1, 0, 1)
init_rho_l <- runif(1, 0, 1)
init_lambda <- runif(1, 0, 10)
init_gamma_gain <- runif(1, 0, 10)
init_gamma_loss <- runif(1, 0, 1)
init_sigma <- runif(1, 0, 1)
solution <- optim(par = c(init_rho_g, 
                          init_rho_l,
                          init_lambda,
                          init_gamma_gain,
                          init_gamma_loss,
                          init_sigma), 
                  fn = cpt_mle, 
                  method = "Nelder-Mead", 
                  control = list("maxit" = 10e6))
print(paste0("Estimated Rho for gains of ", round(abs(solution$par[1]), 2)))

[1] "Estimated Rho for gains of 1.14"

print(paste0("Estimated Rho for losses of ", round(abs(solution$par[2]), 2)))

[1] "Estimated Rho for losses of 0.63"

print(paste0("Estimated Lambda (loss aversion) of ", round(abs(solution$par[3]), 2)))

[1] "Estimated Lambda (loss aversion) of 20.98"

print(paste0("Estimated Gamma for gains of ", round(abs(solution$par[4]), 2)))

[1] "Estimated Gamma for gains of 0.77"

print(paste0("Estimated Gamma for losses of ", round(abs(solution$par[5]), 2)))

[1] "Estimated Gamma for losses of 1"

print(paste0("Estimated Sigma of ", round(abs(solution$par[6]), 2)))

[1] "Estimated Sigma of 30.59"

print(paste0("Convergence check: ", solution$convergence))

[1] "Convergence check: 0"

cpt_predict() predicts Monetary Equivalents for a given dataset:

# - function: CPT model predictions
cpt_predict <- function(data,
                        rho_g, rho_l, lambda,
                        gamma_gain, gamma_loss) {
  
  # - utility functions
  u1 <- rd_utility_power(data$v1, rho_g, rho_l, lambda)
  u2 <- rd_utility_power(data$v2, rho_g, rho_l, lambda)
  
  # - probability weighting
  wp1 <- vector(mode = "numeric", length = length(data$v1))
  ix_gain_1 <- which(data$v1 >= 0)
  wp1[ix_gain_1] <- p_weight(data$p1[ix_gain_1], gamma_gain)
  ix_loss_1 <- which(data$v1 < 0)
  wp1[ix_loss_1] <- p_weight(data$p1[ix_loss_1], gamma_loss)
  wp2 <- 1-wp1
  
  # - predictions: utility scale
  predictions <- wp1*u1 + wp2*u2
  
  # - predictions: value scale
  predictions <- ifelse(predictions >= 0,
                        predictions^(1/rho_g),
                        -(((-predictions)/lambda)^(1/rho_l))
                        )
  
  # - output
  return(predictions)
  
}

NOTE. The predictions: value scale part of the code is very important, because we do not wish to make predictions on the utility scale but rather on the value (MEq) scale. So, let’s not forget that under Prospect Theory we have:

\[MEq = u_g(x)^{1/\rho_g}\] while for losßes we have

\[MEq = -((\frac{u_l(-x)}{\lambda})^{1/\rho_l})\]

Test cpt_predict():

# - predict from the optimized CPT model
meq_data$cpt_meq <- cpt_predict(data = meq_data,
                                rho_g = abs(solution$par[1]),
                                rho_l = abs(solution$par[2]),
                                lambda = abs(solution$par[3]),
                                gamma_gain = abs(solution$par[4]),
                                gamma_loss = abs(solution$par[5])
                                )

# - plot predictions
ggplot(data = meq_data,
       aes(x = meq, y = cpt_meq)) + 
  geom_smooth(method = "lm", size = .25, color = "red", se = FALSE) + 
  geom_point(size = 1.5, color = "black") + 
  geom_point(size = 1, color = "white") + 
  ggtitle("CPT") +
  xlab("Predicted MEq") + 
  ylab("Observed MEq") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

Run 1,000 optimizations (in parallel, using snowfall):

# - init cloud
sfInit(parallel = TRUE, cpus = 23)

R Version:  R version 4.0.3 (2020-10-10) 

# - export
sfExport("meq_data")
sfExport("cpt_mle")
sfExport("cpt_predict")
sfExport("cpt_likelihood")
sfExport("p_weight")
sfExport("rd_utility_power")

# - many optimization runs
init_pars <- lapply(1:1000, function(x) {
  return(
    list(
      init_rho_g = runif(1, 0, 1),
      init_rho_l = runif(1, 0, 1),
      init_lambda = runif(1, 0, 1),
      init_gamma_gain = runif(1, 0, 1),
      init_gamma_loss = runif(1, 0, 1),
      init_sigma = runif(1, 0, 1)
    )
  )
})


# - run
results <- sfClusterApplyLB(init_pars, function(x) {
  
  init_rho_g <- x$init_rho_g
  init_rho_l = x$init_rho_l
  init_lambda = x$init_lambda
  init_gamma_gain = x$init_gamma_gain
  init_gamma_loss = x$init_gamma_loss
  init_sigma <- x$init_sigma
        
  # - optimization run
  optim_sol <- tryCatch(optim(par = c(init_rho_g,
                                      init_rho_l,
                                      init_lambda,
                                      init_gamma_gain,
                                      init_gamma_loss,
                                      init_sigma),
                              fn = cpt_mle,
                              method = "Nelder-Mead",
                              control = list("maxit" = 1e5,
                                             "reltol" = 1e-14)),
                        error = function(condition) {
                          return(NULL)
                        })
  
  if (!is.null(optim_sol)) {
    
    # - predictions
    preds <- cpt_predict(meq_data, 
                         rho_g = abs(optim_sol$par[1]),
                         rho_l = abs(optim_sol$par[2]),
                         lambda = abs(optim_sol$par[3]),
                         gamma_gain = abs(optim_sol$par[4]),
                         gamma_loss = abs(optim_sol$par[5])
                         )
    
    return(
      list(rho_g = abs(optim_sol$par[1]),
           rho_l = abs(optim_sol$par[2]),
           lambda = abs(optim_sol$par[3]),
           gamma_gain = abs(optim_sol$par[4]),
           gamma_loss = abs(optim_sol$par[5]),
           sigma = abs(optim_sol$par[6]),
           nll = optim_sol$value,
           convergence = optim_sol$convergence,
           predictions = preds)
      )
    
  } else {
      
      return(
      list(rho = NULL,
           sigma = NULL,
           nll = NULL,
           convergence = NULL,
           predictions = NULL)
      )
    
  }
  
  
})

# - stop cluster
sfStop()

# - filter out possible degenerative Nelder-Mead solutions
conv_list <- sapply(results, function(x) x$convergence)
degenerate <- which(conv_list != 0)
nll_list <- nll_list[-degenerate]
# - find best solution
nll_list <- sapply(results, function(x) x$nll)
best_run <- which.min(nll_list)[1]
print(paste0("Check convergence: ", results[[best_run]]$convergence))

[1] "Check convergence: 0"

estimated_rho_g <- results[[best_run]]$rho_g
print(paste0("Estimated Rho for gains: ", round(estimated_rho_g,2)))

[1] "Estimated Rho for gains: 0.99"

estimated_rho_l <- results[[best_run]]$rho_l
print(paste0("Estimated Rho for losses: ", round(estimated_rho_l,2)))

[1] "Estimated Rho for losses: 0.8"

estimated_lambda <- results[[best_run]]$lambda
print(paste0("Estimated Lambda (loss aversion): ", round(estimated_lambda,2)))

[1] "Estimated Lambda (loss aversion): 2.36"

estimated_gamma_gain <- results[[best_run]]$gamma_gain
print(paste0("Estimated Gamma for gains: ", round(estimated_gamma_gain,2)))

[1] "Estimated Gamma for gains: 0.6"

estimated_gamma_loss <- results[[best_run]]$gamma_loss
print(paste0("Estimated Gamma for losses: ", round(estimated_gamma_loss,2)))

[1] "Estimated Gamma for losses: 0.77"

estimated_gamma_sigma <- results[[best_run]]$sigma
print(paste0("Estimated Sigma: ", round(estimated_sigma,2)))

[1] "Estimated Sigma: 11.52"

nll <- results[[best_run]]$nll
print(paste0("Negative LogLikelihood: ", round(nll,2)))

[1] "Negative LogLikelihood: 2128.3"

optimal_predictions <- results[[best_run]]$predictions

# - plot predictions
plot_frame <- meq_data
meq_data$cpt_meq <- optimal_predictions
ggplot(data = meq_data,
       aes(x = meq, y = cpt_meq)) + 
  geom_smooth(method = "lm", size = .25, color = "red", se = FALSE) + 
  geom_point(size = 1.5, color = "black") + 
  geom_point(size = 1, color = "white") + 
  ggtitle("CPT") +
  xlab("Observed MEq") + 
  ylab("Predicted MEq") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

PT optimization for the KT92 data set

NOTE. There are no mixed-prospects in the kt92.csv data set:

# - load some MEq experimental data
meq_data <- read.csv(paste0(kt92_data_dir, "kt92.csv"),
                     header = TRUE,
                     check.names = FALSE,
                     row.names = 1,
                     stringsAsFactors = FALSE)
# - percents to probability in meq_data
meq_data$p1 <- meq_data$p1/100
meq_data$p2 <- meq_data$p2/100
print(meq_data)

Run 1,000 optimizations for the KT92 data set:

# - init cloud
sfInit(parallel = TRUE, cpus = 23)

# - export
sfExport("meq_data")
sfExport("cpt_mle")
sfExport("cpt_predict")
sfExport("cpt_likelihood")
sfExport("p_weight")
sfExport("rd_utility_power")

# - many optimization runs
init_pars <- lapply(1:1000, function(x) {
  return(
    list(
      init_rho_g = runif(1, 0, 1),
      init_rho_l = runif(1, 0, 1),
      init_lambda = runif(1, 0, 1),
      init_gamma_gain = runif(1, 0, 1),
      init_gamma_loss = runif(1, 0, 1),
      init_sigma = runif(1, 0, 1)
    )
  )
})


# - run
results <- sfClusterApplyLB(init_pars, function(x) {
  
  init_rho_g <- x$init_rho_g
  init_rho_l = x$init_rho_l
  init_lambda = x$init_lambda
  init_gamma_gain = x$init_gamma_gain
  init_gamma_loss = x$init_gamma_loss
  init_sigma <- x$init_sigma
        
  # - optimization run
  optim_sol <- tryCatch(optim(par = c(init_rho_g,
                                      init_rho_l,
                                      init_lambda,
                                      init_gamma_gain,
                                      init_gamma_loss,
                                      init_sigma),
                              fn = cpt_mle,
                              method = "Nelder-Mead",
                              control = list("maxit" = 1e5,
                                             "reltol" = 1e-14)),
                        error = function(condition) {
                          return(NULL)
                        })
  
  if (!is.null(optim_sol)) {
    
    # - predictions
    preds <- cpt_predict(meq_data, 
                         rho_g = abs(optim_sol$par[1]),
                         rho_l = abs(optim_sol$par[2]),
                         lambda = abs(optim_sol$par[3]),
                         gamma_gain = abs(optim_sol$par[4]),
                         gamma_loss = abs(optim_sol$par[5])
                         )
    
    return(
      list(rho_g = abs(optim_sol$par[1]),
           rho_l = abs(optim_sol$par[2]),
           lambda = abs(optim_sol$par[3]),
           gamma_gain = abs(optim_sol$par[4]),
           gamma_loss = abs(optim_sol$par[5]),
           sigma = abs(optim_sol$par[6]),
           nll = optim_sol$value,
           convergence = optim_sol$convergence,
           predictions = preds)
      )
    
  } else {
      
      return(
      list(rho = NULL,
           sigma = NULL,
           nll = NULL,
           convergence = NULL,
           predictions = NULL)
      )
    
  }
  
  
})

# - stop cluster
sfStop()

# - filter out possible degenerative Nelder-Mead solutions
conv_list <- sapply(results, function(x) x$convergence)
degenerate <- which(conv_list != 0)
nll_list <- nll_list[-degenerate]
# - find best solution
nll_list <- sapply(results, function(x) x$nll)
best_run <- which.min(nll_list)[1]
print(paste0("Check convergence: ", results[[best_run]]$convergence))

[1] "Check convergence: 0"

estimated_rho_g <- results[[best_run]]$rho_g
print(paste0("Estimated Rho for gains: ", round(estimated_rho_g,2)))

[1] "Estimated Rho for gains: 0.62"

estimated_rho_l <- results[[best_run]]$rho_l
print(paste0("Estimated Rho for losses: ", round(estimated_rho_l,2)))

[1] "Estimated Rho for losses: 0.75"

estimated_lambda <- results[[best_run]]$lambda
print(paste0("Estimated Lambda (loss aversion): ", round(estimated_lambda,2)))

[1] "Estimated Lambda (loss aversion): 5.16"

estimated_gamma_gain <- results[[best_run]]$gamma_gain
print(paste0("Estimated Gamma for gains: ", round(estimated_gamma_gain,2)))

[1] "Estimated Gamma for gains: 0.7"

estimated_gamma_loss <- results[[best_run]]$gamma_loss
print(paste0("Estimated Gamma for losses: ", round(estimated_gamma_loss,2)))

[1] "Estimated Gamma for losses: 0.74"

estimated_gamma_sigma <- results[[best_run]]$sigma
print(paste0("Estimated Sigma: ", round(estimated_sigma,2)))

[1] "Estimated Sigma: 11.52"

nll <- results[[best_run]]$nll
print(paste0("Negative LogLikelihood: ", round(nll,2)))

[1] "Negative LogLikelihood: 189.05"

optimal_predictions <- results[[best_run]]$predictions

# - plot predictions
plot_frame <- meq_data
meq_data$cpt_meq <- optimal_predictions
ggplot(data = meq_data,
       aes(x = meq, y = cpt_meq)) + 
  geom_smooth(method = "lm", size = .25, color = "red", se = FALSE) + 
  geom_point(size = 1.5, color = "black") + 
  geom_point(size = 1, color = "white") + 
  ggtitle("CPT") +
  xlab("Observed MEq") + 
  ylab("Predicted MEq") + 
  theme_bw() + 
  theme(panel.border = element_blank()) + 
  theme(plot.title = element_text(hjust=.5, size = 10))

NOTE. \(\lambda\), the Loss Aversion parameter, cannot be estimated from experimental designs that do not encompass mixed lotteries.

Model selection: plan for K-fold CV

# - load some MEq experimental data
meq_data <- read.csv(paste0(meq_data_dir, "es4v1.csv"),
                     header = TRUE,
                     check.names = FALSE,
                     row.names = 1,
                     stringsAsFactors = FALSE)
# - percents to probability in meq_data
meq_data$p1 <- meq_data$p1/100
meq_data$p2 <- meq_data$p2/100
print(meq_data)

Types of lotteries

# - types of lotteries
# - strictly positive
pos_lott <- which(meq_data$v1 > 0 & meq_data$v2 > 0)
print(paste0("There are ", length(pos_lott), " strictly POSITIVE lotteries."))

[1] "There are 90 strictly POSITIVE lotteries."

# - non-negative
nneg_lott <- which ( 
  (meq_data$v1 > 0 & meq_data$v2 == 0) |
    (meq_data$v1 == 0 & meq_data$v2 > 0)
)
print(paste0("There are ", length(nneg_lott), " NON-NEGATIVE lotteries."))

[1] "There are 45 NON-NEGATIVE lotteries."

# - strictly negative
neg_lott <- which(meq_data$v1 < 0 & meq_data$v2 < 0)
print(paste0("There are ", length(neg_lott), " strictly NEGATIVE lotteries."))

[1] "There are 90 strictly NEGATIVE lotteries."

# - non-positive
npos_lott <- which ( 
  (meq_data$v1 < 0 & meq_data$v2 == 0) |
    (meq_data$v1 == 0 & meq_data$v2 < 0)
)
print(paste0("There are ", length(npos_lott), " NON-POSITIVE lotteries."))

[1] "There are 45 NON-POSITIVE lotteries."

# - mixed
mixed_lott <- which ( 
  (meq_data$v1 < 0 & meq_data$v2 > 0) |
    (meq_data$v1 > 0 & meq_data$v2 < 0)
)
print(paste0("There are ", length(mixed_lott), " MIXED lotteries."))

[1] "There are 225 MIXED lotteries."

Let’s enter type to meq_data:

meq_data$type <- apply(meq_data, 1, function(x) {
  if (x[2] > 0 & x[4] > 0) {
    return("p")
  } else if ((x[2] > 0 & x[4] == 0) | 
             (x[2] == 0 & x[4] > 0)) {
    return("nn")
  } else if (x[2] < 0 & x[4] < 0) {
    return("n")
  } else if ((x[2] < 0 & x[4] == 0) |
    (x[2] == 0 & x[4] < 0)) {
    return("np")
  } else {
    return("mixed")
  }
})
table(meq_data$type)

mixed     n    nn    np     p 
  225    90    45    45    90 

You need to make sure that all types of lotteries enter all of your folds in k-fold cross-validation of this model.

Folds plan

meq_data$fold <- 0
types <- unique(meq_data$type)
for (i in 1:length(types)) {
  w <- which(meq_data$type == types[i])
  f_asgn <- rep(1:5, length(w)/5)
  f_asgn <- sample(f_asgn, size = length(f_asgn), replace = FALSE)
  meq_data$fold[w] <- f_asgn
}
print(meq_data)

Check:

table(meq_data$fold, meq_data$type)

   
    mixed  n nn np  p
  1    45 18  9  9 18
  2    45 18  9  9 18
  3    45 18  9  9 18
  4    45 18  9  9 18
  5    45 18  9  9 18

---

Goran S. Milovanović

DataKolektiv, 2022.

contact: goran.milovanovic@datakolektiv.com

---

License: GPLv3 This Notebook is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This Notebook is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this Notebook. If not, see http://www.gnu.org/licenses/.

---